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1 Introduction 

Description of simple hadrons in terms of quark-gluon degrees of freedom has long been 
an active area in physics. With the advent of TJNAF, which will be operating at inter- 
mediate energies and therefore probing the structure of hadrons, there is new motivation 
and need for a simple theoretical description of quark interactions. In this context, the 
Dyson-Schwinger Bethe-Salpeter(DSBS) equation formalism has gained popularity in re- 
cent years.Q The DSBS formalism serves to bridge the gap between nonrelativistic quark 
models and more rigorous approaches, such as lattice gauge theory. 

The main features of QCD can be summarized as chiral symmetry breaking, confine- 
ment and asymptotic freedom. It is possible to address all of these features within the 
DSBS formalism. In this formalism, the input is an effective gluon propagator which is 
assumed to represent the interactions between quarks at all momentum transfers. The 
choice of a vector interaction between quarks is motivated only by the desire to make a 
connection with QCD degrees of freedom. In fact, whether a scalar or a vector interaction 
should be used between quarks is a topic of debate not addressed in this paper. 

While various applications of the DSBS formalism to pseudoscalar and vector mesons 
have produced promising results, there are still some questions to be investigated. In this 
paper, we address three issues. These are: a) Can an infrared gluon propagator lead 
to confined quarks? b) The question of using Euclidean metric and extrapolation, c) 
Dressing of the quark-gluon vertex in the DSBS equations while maintaining the chiral 
limit ? 

The organization of this paper is as follows: In section H, the model is introduced 
and the realization of the chiral limit, within the dressing scheme used in this paper, is 
discussed. In section ^ the quark propagator functions obtained by solving the Dyson- 
Schwinger equation are presented and the quark propagator is shown to be free of real 
timelike poles, indicating that quarks can not be free, which is an implication and require- 
ment of confinement. In section ^, meson and kaon observables, which are calculated using 
a Euclidean metric(rather than a Wick rotation that only effects the internal momenta), 
are presented. Finally, results are summarized and our conclusions are presented in section 



■^See Ref. ^ for an extensive review. 
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Figure 1: Quark Dyson- Schwinger equation is shown. Only one vertex is dressed to 
prevent double counting. 
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2 THE MODEL 

The mechanism of chiral symmetry breaking and recovery of massless pseudoscalar bound 
states(pion, kaon) in the limit of massless fermions( quarks) was originally discovered in 
the papers of Nambu-Jona-Lasinio(NJL) [^. Nambu-Jona-Lasinio's model originally de- 
scribed relativistic nucleon interactions through local, four-nucleon couplings. It is the 
same philosophy that is followed in the DSBS calculations, except that now nucleons 
are replaced with quarks and the contact interaction is replaced by an effective gluon 
exchange between quarks. The Dyson- Schwinger (DS) and the Bethe-Salpeter(BS) equa- 
tions employed in this work are shown in Figures |I| and ^. The DS equation describes 
the propagation of quarks in the presence of gluons. The BS equation describes the 
quark- ant iquark bound state, in which the DS quark propagator is used. 

In Fig. in, the thin lines represent the current quark propagators for each flavor f: 



^^?(p)=[-7 ol (1) 

while the thick lines correspond to the dressed quark propagators: 



Here Af{p) is a dimensionless normalization factor and Bf{p) has the units of mass(MeV). 

The problem of how to systematically dress DS and BS equations has recently been 

addressed[^, |^. Dressing of all vertices consistently is motivated by the desire to make 

a closer connection with QCD. The structure of the Dyson-Schwinger equation, when 




Figure 2: The Bethe-Salpeter equation is shown. Only one vertex is dressed each time. 
This is necessary to preserve the chiral hmit, which follows from the similarity of the BS 
and the DS equations. 

combined with the chiral limit requirement, strictly restricts the choice of kernel for the 
Bethe-Salpeter equation. The chiral limit in NJL type models f\ is obtained due to the 
similarity of the Dyson-Schwinger and Bethe-Salpeter equations in the limit of massless 
current quarks. In this limit, the quark mass function B{p) and the Bethe-Salpeter wave- 
function $(p) for pseudoscalar massless bound states satisfy the same equation(to be 
shown below). Therefore, for any given set of parameters of the model gluon propaga- 
tor G^'^ , the solution of the Dyson-Schwinger equation automatically implies a massless 
pseudoscalar bound state solution for the Bethe-Salpeter equation. In other words, in the 
chiral limit the quark Dyson-Schwinger equation produces the appropriate mass function 
B{p) such that the BS equation produces a massless pseudoscalar bound state. 

In order to prevent double counting, only one of the vertices in the Dyson-Schwinger 
equation(Fig. ||) is dressed as indicated by the solid circle. Therefore, to preserve the 
similarity between the BS and the DS equations in the chiral limit, we dress only one of 
the quark-gluon vertices in constructing the BS equation. In order to keep the Bethe- 
Salpeter equation symmetric(to treat quarks equally), the kernel is divided into two pieces, 
where in each piece an alternate vertex is dressed, and contribution of those terms is 
averaged(See Fig. ^. While the dressing of only one of the quark-gluon vertices in the BS 
kernel does not represent a complete dressing, since cases where both quark-gluon vertices 
are simultaneously dressed are excluded, with the proper choice of vertex F'^, one has a 
subset of all diagrams that produce the correct chiral limit. 

Having stated the general structure of the DS and BS equations used in this calcula- 

*See Refs. p, n for an extensive review of the NJL type models, and Refs. [bI for an extended version 
of it. 



tion, we now discuss the chiral limit, and the choice of the gluon propagator G^^{q). 
In terms of quark and gluon propagators, the DS equation is written as 

S}{p) = Sofip) + ^lSf{p)[J ^ G^-'ip - q) r,(p, q)Sf{q)^AS,j{p). (3) 

Similarly, the Bethe-Salpeter equation determining the BS vertex function(a truncated 
wavefunction, see later) Xp(^) is given by 
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where the 4- vector q^ = Prji + q, q- = Pr]2 — q, rji + ri2 = 1, and P is the bound state 
4- momentum. The BS vertex function Xp{k) and its conjugate Xp(^) are related [l^ by 



Xpiiko, k) = 7oXp(-^^o, fc)7o- 



(5) 



The normalization condition of the BS vertex function is derivable from the BS equation 
itself. [|11|] If the interaction does not depend on the total bound state momentum, the 
normalization condition reduces to 
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As it is well known, the dressing of electromagnetic vertices, such as the photon-quark 
vertex, is constrained by the Ward-Takahashi identity. 



q,T'^{p',p) = S-\p')-S-\p), 



(7) 



which guarantees the conservation of electromagnetic current at the vertex. Similarly, 
due to color current conservation, the dressed quark-gluon interaction vertex, r'^(p, g), 
satisfies the Slavnov- Taylor identity [14| 



q,T^{p',p)[l + b{q')] = [1 - Biq,p)]S-\p') - S-\p)[l - B{q,p)], 



where functions &(g^) and B{q,p) are related to "ghost fields." Since the QCD ghost field 
contributions to Eq. |^ are not well understood we neglect their contributions. When ghost 
fields are neglected, 6(g^) and B{q,p) = 0; the Slavnov- Taylor identity then reduces to 
the Ward-Takahashi identity Eq. |^. The minimal vertex that satisfies this identity in this 
limit has been given by Ball-Chiu [|15] as: 
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It is clear that one can add any term to this vertex that satisfies 

q,rf^{p',p) = 0. 



(f+^) + Bip')-Bip) (9) 



(10) 



Curtis-Pennington |T6| have proposed such an additional vertex term. Here, for simplicity, 
we consider only the Ball-Chiu dressing. 



2.1 The Chiral Limit 

To determine what type of quark-gluon vertex is allowed within the approximation scheme 
employed here, let us analyze the chiral limit of the DS and the BS equations. In the 
chiral limit the current quark masses ttlq and the bound state momentum P vanishes. In 
this limit, taking the spinor trace of the Dyson-Schwinger equation |^, one obtains 
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For a pseudoscalar bound state, in the chiral limit the BS vertex function is given by 



Xp(fc) = Z75*p(fc), 



(12) 



where ^p{k) is a scalar function of the relative momentum k, and bound state momentum 
P = 0. With this definition, the BS equation for the bound state wavefunction $p(A;) is 
obtained as 

i r d^q 
3 
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l,T''{k,q)SMl^^p{q)Su{q)l, (13) 
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Commuting 75 through propagators and noting that in the chiral hmit, 

Eq. |l^ can be rewritten as 

^,ik) = '-( ^^G^^iq-k) ^-^ tr[(r,(g,fc)-75r,(fc,g)75)% 

To ensure that B{k) and ^p{k) satisfy the same equation in the chiral hmit, therefore 
producing a massless pion, the following condition should be satisfied 

..[r,.fe..„i.i|„-,„] . "['^-'^■^';^-"'-^»-'-' ,15) 

Clearly, the commonly used bare quark gluon vertex, F^ = 7^ satisfies this condition. 
Furthermore, the condition Eq. |T5| is also satisfied by the following part of the Ball-Chiu 
vertex 

r.(p'.rt . m±m.,^ii±iiim^^^-,^y ,16) 

The -B(p') — B{p) part of the Ball-Chiu vertex would have contributed to the lefthandside 
of Eq. |T3] while it does not contribute to the righthandside. The contribution of this term 
is small compared to the sum of other two terms. Therefore, it is a good approximation to 



drop this term. The vertex Eq. |T6| enables us to preserve the chiral limit while including 
the dominant part of the quark-gluon vertex dressing as suggested by the Ward-Takahashi 
identity. This approximate form(Eq. |1^) of the Ball-Chiu vertex is used for the quark- 
gluon coupling in this work. 

Since the functions B{k) and ^p{k) satisfy the same equation in the chiral limit, 
they are equal upto a proportionality constant. With the help of the BS normalization 
condition |^ it is found 0, |l^ that 

^p{k) = ^B{k), (17) 

where /^ is the pion decay constant. 

In order to complete the description of the model, one needs to choose a model for the 



gluon propagator G^''{q). The choice of G^''{q) has been discussed in various papers jT? 



Usually, the ultraviolet (g -^ oo) or asymptotic behavior of G^^{q) is borrowed from 
QCD calculations, while its infrared(g — > 0) or confining behavior is given by a sharply 
falling function such as 1/g^ or ^^(i?), to incorporate confinement. The problem with 1/q^ 
behavior is that it is not an integrable singularity and one needs to introduce an infrared 
cutoff. While the 5^{q) form does not have this problem, it is perhaps too simple a form 
to represent the physics in the infrared region; this issue requires further study. In this 
paper, we show that an infrared finite propagator can not be ruled out on the 
basis of quark Dyson-Schwinger and Bethe-Salpeter equations. The model used 
here is 

G^'iq) = ig^" - ^) [GiR + Guv], (18) 

where the infrared(g -^ 0) behavior, Gjr, is modeled by a finite term 

GiR(g) = Ge-«'/-', (19) 

while the asymptotic(g -^ oo) form. Guv, is taken from perturbative QCD calculations 

with d = 12/(33 — ^Nf) = |, where Nj = 3 is the number of flavors. The QCD scale 
parameter Aqcd, determined by fitting high energy experiments(Particle Data Group, 
1990), is chosen to be 225 MeV. The constant r ensures the positivity of the asymptotic 
piece as g ^ 0, and results are not very sensitive to this parameter. The r is chosen to 
be r = 3. The a^ in the denominator of the UV piece is introduced to ensure that the 
infrared piece is the dominant contribution at low energies. Similar forms for the dressed 
gluon propagator have been used in the literature |]TH|, [T^, [2^, ^ with considerable success 
in preliminary applications of Dyson-Schwinger Equations to hadronic physics. Here, we 
develop this approach using a dressed quark-gluon vertex while maintaining the chiral 
limit, and show that quarks can be confined with an infrared finite interaction. 



Aside from current quark masses, m^.d ~ 6 ± 2 MeV and m^ ~ 150 ± 50 MeV |22[, which 
are also the input parameters of QCD, there are only two unconstrained parameters(G, cr) 
to vary to predict the data. Parameters G = 1.9710~^MeV~^ and a = 750 MeV are chosen 
to give the optimum overall fit. We choose the current quark mass of the strange quark to 



fine tune the kaon mass. The current quark masses used in the calculation are: rriu^d = 3 
MeV, and rus = 60 MeV. The ratio nis/niu is well within acceptable limits ||22|. In the 
next section, we discuss the solution of the quark Dyson-Schwinger equation. 



3 QUARK PROPAGATORS AND CONFINEMENT 

The numerical solution of the DS equation Eq. ^ is performed through iteration to find 
the quark propagator functions A{p'^) and BijP) in Euclidean metric. The details of the 
numerical methods are explained in the Appendix. Solutions for A{p'^) and M{p'^) = 
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Figure 3: The quark propagator function A{p) is shown for up/down, and strange quarks. 
Dressing of the quark-gluon vertex causes the peak observed in A{p). 

B{p^) / A{p'^) in the spacelike region are shown in Figures |]and|^. According to the results 
in Figures |^ and ^, for low momenta (p < 1 GeV) quark masses are close to those of the 
constituent quarks; whereas, as momentum increases, quarks start behaving as current 
quarks(y4/(g^) — > l,_B/(g^) — * mo). In the infrared region, the quark masses approach 
the constituent mass values. The behavior of function Af{q^) is different than that of 
Mfijp)] Af{q^) reaches its maximum value at intermediate energies, where the scale is 



given by a. This is a different behavior than results presented in Ref. [0, for the bare 
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Figure 4: Quark mass functions M{p) are shown for up/down, and strange quarks. At 
the origin(p^ = 0), quark masses are closer to constituent quark mass values used in 
nonrelativistic quark models. Asymptotically, quark mass values approach current quark 
mass values. 

quark-gluon coupling case. In their results, both M(p^) and Aj(g^) are monotonically 
decreasing functions. In our case, the dressing of the quark-gluon vertex gives rise to the 
nonmonotonic behavior we found for Af{q'^). As one increases the coupling strength G, 
monotonic behavior of A is restored. 



3.1 Test of Confinement 

Confinement is the property that only color singlet hadrons are observed in nature. In 
QCD, confinement is obtained dynamically due to the nonabelian, hence self-interacting, 
nature of gluons. A natural result of confinement is that no free quark state should be 
observed(a free quark state has a net color.) In QFT, an n-body bound state is defined 
by the pole of the n-body propagator. The familiar 2-body(Bethe-Salpeter, Gross) and 
3-body(Faddeev) equations are obtained, based on this definition, by looking for the poles 
in the two and three body propagators. Similarly, it is natural to expect that a l-body(or 
free) state should be identified by the pole of the one body propagator. Therefore, if 



quarks are confined, quark propagators should not have poles in the tinielike(p^ > 0) 
region.0, because such a pole permits an asymptotically (in space) free quark wave to 
exist. The absence of poles in the quark propagators is, however, not conclusive evidence 
for confinement. In order to be able to claim that a theory is confining, one has to also 
show that a diquark, or any other color nonsinglet stable bound state does not exist. 
Here, we restrict our discussion of confinement to quarks only. Therefore from here on, 
"confinement" refers to "lack of free quarks" rather than the more general definition of 
"lack of colored states." In order to test whether a quark propagator, given in Euclidean 
metric, leads to confinement, one needs a procedure to determine the presence of any poles 
in the timelike region.^ For simple cases such as a free fermion, the Euclidean expression 
for the propagator can be readily used to see if there are any poles when the propagator 
is continued to the Minkowski metric, i.e. 

where there is a pole at p^ = m^. For the dressed quark propagator, the confinement test 
question is: does 

^(P^' (22) 



^'(pIH + BHA) 



have a pole? For this test, the procedure used in Ref. |2^ is adopted to determine whether 
a quark propagator given in Euclidean metric has poles, when continued to Minkowski 
metric. The starting point is the definition of a generalized quark propagator 

S{p)^ fdf,p{f,)-^, (23) 

J p- fi 

where p(/x) is a spectral density function. The Euchdean metric expression for tr[S'£;(p)] 

is 

/i 



tr[5^(p)]=4 dfipif,)^^. (24) 

p"^ + jj,^ 



''There is an alternative to this approach, which is developed within the context of the Gross equation 



in Ref. |23 . In that approach, quarks are allowed to be on shell as long as they are in the vicinity of 
off-shell quarks, and confinement is realized through a relativistic generalization of the linear potential. 
^An alternative realization of confinement can be obtained by simply defining the quark mass function 
such that quark can never be on shell. g This approach amounts to having quark mass function Af(p^) 
as input and the effective gluon propagator G^'^{q) as unknown in the quark Dyson-Schwinger equation. 
In this approach, a unique determination of the gluon propagator is not possible. Therefore, one is forced 
to make a separable interaction approximation. 

10 



Using contour integration methods, the Fourier transform of this trace is found to be 

A(t) = J^e^'hrlSEip)] (25) 



2 f dfx pifx) e-f"^ . (26) 



We define /i(t), the average effective mass(pole location) of the propagator S{p) as a 
function of time by 

m = -|i^(^w) (27) 

J dfi p{fi) e-t^' 
Let us assume that S{p) has at least one pole for a finite p"^ > 0. Let p{fi) = J2i=o Cj S{p — 
rrii), where n > 1 and rrii < rrii^i. This assumption leads to a discrete average 






M^) = ^li\ '. • (29) 



Taking the limit of this expression at infinite time t ^ oo, one has 

lim p{t) ^ mo. (30) 

t — >oo 

Therefore, this averaging procedure gives the smallest pole of the propagator S{p). If this 
limit exists and it is real, then there is at least one finite pole for timelike momenta in the 
Minkowski metric, and therefore the propagator does not represent a spatially confined 
particle. On the other hand, if there is no finite limit or the limit is complex, then the 
propagator represents a confined particle, since the propagator does not have a real finite 
pole. This test for confinement has been applied to the quark propagator obtained from 
the numerical solution of the Dyson-Schwinger equation. There are three possible pole 
structures for the quark propagator. These possibilities are exemplified by the following: 

• Real pole: 

1 



SsiPi 



p^ + m'^ 
A(t) oc e"'"*, 

pit) = m, 
where the test produces the pole location as expected. 
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Absence of poles: 



A(t) oc e^''"'"l\ 
p,{t) = crt — > oo, 

In this example, there is no finite pole and the propagator represents a confined particle. 
This analytic form is the same as the infrared piece of the effective gluon propagator |T9|. 

• Complex poles: 

1 



Se{pi 



(Pe ~ io,)'^ — Tn-^ 

Here the pole is complex in general and purely imaginary when a = 0. For this case, the 
test results gives analytically 

A(t) oc e""*sin(mt), 

/2(t) = a + mtan(mt — — ). 

Therefore, the signature of complex poles appears as M(t) = ]i(i) oc tan(mt) behavior as 
t ^ oo, where the frequency of oscillations is proportional to the imaginary part of the 
quark mass pole. This is exactly the type of behavior found by applying the confinement 
test to the quark propagator obtained by numerically solving the DS equation ^. Since 
the quark propagator is known only numerically, application of the test is numerical 
and details of the numerical methods are explained in the Appendix. Here we present 
the test results for three different coupling strengths, namely G = 1 x 10~^, 1.5 x 10~^, 
and 1.97 x IQ-^MeVl The first case, G = 1 x IQ-^MeV"^ is shown in Figure |. 
According to this result, the pole location is finite(~ llOMeV) and real. Therefore, 
this quark propagator does not represent a confined particle. On the other hand, this is 
a case where the coupling constant is very small. As one increases the strength of the 
coupling to G = 1.5 10~^MeV~^ an irregular oscillatory behavior sets in(Fig. ^. If the 
coupling strength is further increased to G = 1.97 10~^MeV~^, which is the parameter 
that is used to fit all observables in this paper, the oscillations clearly displays the 
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Figure 5: Mass pole as a function of time is shown. Asymptotically, M(t) approaches a 
constant, which indicates a real mass pole(unconfined quark). 



tan(mt) behavior(Fig. 0) which indicates that quarks have complex mass poles. 
According to this result (Fig. |^, the average distance a quark can travel before it 

hadronizes, which is given by the average distance between the peaks (singularities) of the 
M{t) function, is approximately D x 200MeV/cr Fermi = 1.45 x 200/750 = 0.39 Fermi, 
where D ^ 1.45 is the average spacing between the peaks. This result is in very good 
agreement when compared with the sizes of various hadronic bound states such as the 
pion(r,r = -66 Fermi) and the kaon(rx = -53 Fermi). A recent study |2^ of the DS 
equation within the context of QED in three dimensions similarly finds complex conjugate 
mass poles in the fermion propagator. It is important to emphasize that the analytic 
form of the gluon propagator used in this calculation is not the only possible choice to 
produce confined quarks. In fact, we have obtained similar results with other infrared 
singular gluon propagators. Therefore, it is not possible to single out a specific analytic 
form for the gluon propagator solely on the basis of the confinement test. Just as there 
are infinitely many confining potentials in nonrelativistic quantum mechanics, there are 
infinitely many confining effective gluon propagators in field theory. Therefore, additional 
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Figure 6: Mass pole as a function of time is shown. Oscillations indicate that quark mass 
pole is complex. 
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Figure 7: As the strength of the coupling, G, is increased M{t) oc tan(t) behavior clearly 
sets in. 
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tests such as the prediction of meson observables are needed to determine if the gluon 
propagator ansatz makes physical sense. 

Having shown that it is possible to obtain confined quarks using a gluon propagator 
with a gaussian type of infrared behavior, we now turn to the quark-antiquark bound 
state problem. 

4 QUARK-ANTIQUARK BOUND STATES 

Before we embark on solving the BS equation, it is necessary to clarify a technical prob- 
lem. In the previous section, dressed quark propagators were calculated in the Euclidean 
metric(or spacelike momentum region). For the Bethe-Salpeter equation, usage of the 
Euclidean metric is more problematic .[] Unlike the quark Dyson-Schwinger equation, the 
Bethe-Salpeter equation involves the external total bound state momentum P, which has 
to eventually represent a physical particle(bound state) with a real positive mass. There- 
fore, the four momentum of the particle should be P = (m, 0), for which P^ = m? > 
represents a timelike particle. It follows that, in order to be able to perform the inte- 
grations in Eq. ^ in Euclidean metric, one needs to know the functions A{q'^),B{q'^) for 
q1_ = rn?rjirj2 — q^ + imqo, where 771^2 is chosen to be 

^1,2 = f , 31 

rrii + 1712 

and mi^2 = Mu,d{0). On the other hand, functions yl(g^), i?(g^) are known only for real 
and spacelike g^. The thick line in Fig. ^(the positive g^ axis) is where functions A{q'^) 
and B{q^) has been calculated by solving the Dyson-Schwinger equation. The domain 
where these functions are needed is shown by the shaded region in Fig. p. At this point, 
there are three options. The first one, used in Ref. |^ among other works, is to assume 



an analytic functional form that fits the numerical functions A(g^), B{q^) on the positive 
real axis. Once this assumption is made one can use these analytic functions over the 



entire complex plane. The second approach, which has been used in Ref. |]T8[, is to Taylor 
expand the functions A{q^),B{q^) around their real values to extend the results to the 
complex plane. Both of these methods have drawbacks. It has been shown in previous 
^Recent efforts to solve the BS equation direetly in Minkowski metric are given in Rcfs. pO, M 
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Figure 8: Argument domain of functions A(q'+) and B{q^) for physical bound states(r7i^ > 
0). The positive axis shows the spacehke q\ region. 

section that the quark propagator has complex conjugate mass poles in the spacelike 
region. It has also been shown p9[| that functions A{q^)^B{q^) possibly have poles and 
branch cuts in the complex plane, which complicates the above methods. Therefore 
sampling the momentum values away from the spacelike momentum axis is problematic. 
We follow a third approach, which is also used in lattice gauge theory calculations. In this 
approach, the bound state problem is solved for a set of spacelike bound state momenta, 
P^ = —rn^ < 0, which transforms the problem to Euclidean metric, thereby avoiding 
the complex argument problem, and the final results are then extrapolated back to the 
physical region, P^ > 0. Since it is the final results such as form factors and decay 
constants that are extrapolated ( rather than functions A{q'^), B{q^), which are integrated 
out in calculations of observables) , this method has the benefit of explicitly displaying the 
reliability of the extrapolation. The only assumption is the analyticity of the observables 
as a function of the bound state mass. In this method, there is no need to assume 
that functions A(g^),i?(g^) are analytic. The procedure for the solution is as follows: 
First the Bethe-Salpeter equation ^ is discretizedQ and transformed into a matrix Equation 



Hp2 $ = $, 



(32) 



where P^ = rn^ > is the bound state mass and the implicit eigenvalue of this matrix 
equation. Since Eq. ^ will be solved for spacelike bound states P^ = ml < 0, one will 
^see Appendix for details of our numerical techniques. 
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not be able to find any solutions unless an artificial eigenvalue ai is introduced to Eq. p2| 

H^2 $ = ai^. (33) 

One proceeds by finding a set of solutions {mf < 0, aj} to the above equation. This is 
done by using an inverse iteration technique as explained in the Appendix. Using the set 
of solutions {m"^ < 0, Oj}, a functional relationship between ctj and mf can be established 

a, = /(m^). (34) 

It is only when mf = m? > and a^ = 1 that one recovers the original BS equation(^2|). 
Therefore, location of the rn^ that gives a = 1 is the eigenvalue and mass of the physical 
bound state in question. The most general form for the spin-space part of the BS vertex 
function for pseudoscalar mesons is given by 

Xp(fc) = 275[$0 + f $1 + ^ $2 + [ 1^, f ] $3]. (35) 

In the pseudoscalar meson channel, the dominant contribution to the BS vertex function 



comes from the first term [0, [T^ 



Xpik) ^ «75$p(A;) (36) 

This dominance is not surprising since the Dirac structure of the leading term in ^ is the 
same as that of pointlike pion-quark coupling. Therefore, we only consider the leading 
term in our analysis.^] The angular dependence of the BS vertex function is made explicit 
by expanding it in terms of Tchebyshev polynomials 

oo 

$p(A;) = 5:$^(|fc|)ncos7). (37) 

n=0 

For the bound states considered in this work(7r, vr*, K^), the dominant contribution to 
$p(A;) comes from the Tq polynomial. Test runs for cases where higher(n > 0) Tcheby- 
shev polynomials are included showed that the contribution of the higher Tchebyshev 
polynomials are negligible, which is in agreement with conclusions in Ref. |T^. After 



'For vector mesons, the leading term might not be the only important one. 
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Figure 9: Dependence of eigenvalue a on P^ 
momentum is given by the intersection defined by a 



m^. Location of the bound state 



solving the BS equation numerically, a relationship |^ between the largest eigenvalue^ ot 
and vp? is constructed by using a polynomial fit 



j=0 



(3^ 



It has been determined that for all of the bound states under consideration n = 3 gave 
a satisfactory fit. The extrapolations done to find the pion and kaon masses are shown 
respectively in Fig. and Fig. [TO. The pion ground and first excited state and kaon 



groundstate wavefunctions are shown in Fig. [T^. The first excited state of the pion 
has a node, while the ground state wavefunctions are positive definite. The kaon BS 
wavefunction is more spread out in momentum space than that of pion. 

4.1 The Electromagnetic form factor using Euclidean metric 

Electromagnetic current conservation 



(9 7^ = n 



(39) 



10 



The second largest a leads to the first excited state. 
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Figure 10: Determination of kaon mass. 
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Figure 11: The pion and kaon ground and pion first excited state BS wavefunctions 
$p(|/i;|) are shown all for P^ = 0. 
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photon 




Figure 12: Electromagnetic Form Factor of Pion in the impulse approximation(i.e. inter- 
action of the initial and final pion with each other is neglected.) 

implies that in momentum space one must have Q'/^J^^ — 0- Since both initial and final 



pion states are on their mass shells, we have p"^ 



m 



ij) + g)^, or 2p ■ q + q^ 



0. 



Therefore, the definition of the form factor takes the following form 



<7r+(p + g)|J,^JO)| 7r+(p) >= FM'){'^P' + Q') 



(40) 



This matrix element is represented by the Feynman diagram shown in Figure |T2[^ As 
the photon momentum q vanishes, the pion is perceived as a point charge. Therefore, it 
is expected that -F(O) = e, and the charged pion has an electromagnetic charge of e. To 
arrive at this result, two conditions must be satisfied; namely, the BS wavefunction should 
be properly normalized, and the conservation of the electromagnetic current at the photon- 
quark interaction vertex should be taken into account. The dressed electromagnetic vertex 
function r^(p',p) for the photon quark coupling satisfies the Ward-Takahashi identity 
Eq. 1^, which is an indirect statement of the conservation of the electromagnetic current. 
For the photon-quark vertex we use the full Ball-Chiu vertex ^ In the Breit frame, the 
initial and final pion momenta are 



Pi 



(^v/"^?- j,0,-|), P/ = (iWm2-^,0,|). 



(41) 



Therefore, (-Pj)^ = (-P/)^ = —fn^ < are both spacelike. Calculation of the form factor is 
done for a set of spacelike bound state masses mf, {i = 1 ■ ■ -21) and the resultant set of 
form factor functions are extrapolated at each photon momentum. For example, the pion 



^^Only the interaction of the photon with the u quark is shown. There is also a second diagram where 
the photon interacts with the d quark. 
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Figure 13: These plots show how the form factor extrapolation is done for various photon 
energies. The extrapolated results are on the right hand side(tinielike, P^ = m^ > 0) of 
bound state momentum axis. 



form factor at photon momentum of qf is 



j=0 



2\ 27 



(42) 



where C^{qf) are the coefficients of a fit to the numerical result for F for photon momen- 
tum of qf. Therefore, the physical result for the pion form factor at the photon momentum 
qf is F{M^,qf). In Figure |l^, we show the extrapolation for three different photon mo- 
menta, g2 = 0.085 GeV2, q^ = 0.16 GeV^, and g^ = 0.275 GeV^. For each case, the result 
of extrapolation is given by the value of the form factor functions at P^ = rn^. In order 
to ensure the reliability of the fit, a large number(21) of form factor calculations at each 
photon momentum has been performed. The order of the polynomial fit used is n = 7. 
According to the extrapolation results(Fig. [13|) as the photon momentum increases, the 
reliability of the extrapolation decreases. This difficulty is common to all solution meth- 
ods(within the covariant DS-BS formalism) which use Euclidean metric.0 The reason for 

""^^See Refs. p3] for a discussion of the transition to asyinptotic(large q^) behavior of the pion form 
factor within a hght-cone Bethe-Salpeter forniahsm. 
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Figure 14: Extrapolated result for electromagnetic pion form factor. 



this sensitivity is the increase in the curvature of the function F{m'^,qf) with increasing 
qf. We present the form factor calculation for pion(Fig. |14D and kaon(Fig. |15D cases in 
the region where the extrapolation is reliable. The data points for the pion form factor 
are taken from Refs. ^^ |3^. Data available from earlier experiments [^, ^ for the kaon 
form factor are very poor. Experiments at TJNAF will hopefully provide better measure- 
ments for both pion and kaon form factors. The pion decay constant, /^r, is defined by 
the vacuum to one pion matrix element of the axial vector current: 



<0\Jl,{x)\7i^{p)>^tfJ'^p,e-'^\ 



For a TT^ meson at a; = 0, this definition becomes: 



(43) 



< 0|fr,(0)7^75 — ^,(0)|7r+(p) >= tf^p^. 



(44) 



This matrix element corresponds to the Feynman diagram shown in Figure |T^. Calculation 
of decay constants are done using the same type of extrapolation technique used in the 
calculations of masses and form factors. Values found are fx = 90.3(113) MeV, and /t^ = 
77.2(92.4) MeV, where numbers in parenthesis are the experimental measurements. The 
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Figure 15: Electromagnetic form factor of Kaon. 
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Figure 16: Pion Decay 
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Figure 17: Neutral Pion Decay. 
ratio of decay constants JkI j-n = 1.17(1.22) is in good agreement with the experimental 



value and comparable to those found in similar works. It is reported in Ref. |^ that 
next to leading order terms of the vertex function |3^ contribute significantly to the decay 
constants. This is a possible explanation for the low values we find for the decay constants. 
As a final application, the neutral pion decay to two photons, tt ^ 7 + 7 is considered. 
Neutral pion decay is of historical importance since it is associated with the axial anomaly. 
The matrix element for vr'' -^ 77 decay(Fig. |l^ is 



T{ki,k2) = -2^3:pe^i.p^ 



xe^*{h) 6''*{k2) kip k2. M{ki, k2), (45) 



Since final photons are on-shell, kf = 0, and P^ = {ki + /c2)^ = 2 ki ■ k2 = m^. Therefore 
the scalar function M{ki, ^2) is only a function of the pion mass. The decay rate 

Ko^,, = i^?^M{ml), (46) 

is experimentally measured as F^o^^^ = 7.74 ± 0.56 eV, which means 

g^o^^ = M{ml) = 0.504 ± 0.019. (47) 

It has been shown in Ref. [^ that, in chiral limit, irrespective of the details of quark 



propagators, as long as the photon-quark vertices are properly dressed to conserve the 
electromagnetic current and the BS vertex function is properly normalized, gj^o-y-y is an- 
al5d;ically found to be 0.5. When the mass of the pion is taken into account, we find 
g^^o^^ = .43 which is close to the experimental value. 

A summary of the observables we calculated is given in Table |I|. Error bars in exper- 
imental measurements are negligible unless indicated. 
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Observable 


Calculated 


Experimental 


m^( MeV) 


148 


139.6 


M MeV) 


77.2 


92.4 


< rl >i/2( Fermi) 


.65 


.66 


gTvO-yy 


.43 


.504 


mli MeV) 


1245 


1300± 100 


mxi MeV) 


515 


495 


/i.(MeV) 


90.3 


113 


< rl >i/2( Fermi) 


.54 


.53 


Ik/ fn 


1.17 


1.22 


rK/r\ 


0.83 


0.8 



Table 1: Summary of results 

5 CONCLUSION 

In this work, we analyzed three aspects of the Dyson-Schwinger Bethe-Salpeter equation 
approach. We have shown that it is possible to dress the quark-gluon vertex beyond 
the simple 7^^ form while maintaining the chiral limit condition. We have shown that 
an infrared finite gluon propagator can lead to confined quarks through generation of 
complex quark masses. It was found that, according to the model presented here, quarks 
can freely propagate only ^0.4 Fermi which is in very good agreement with the hadronic 
bound state sizes. We have calculated all observables, including the pion form factor, 
using a Euclidean metric approach without relying on the analyticity properties of quark 
propagator functions A{p'^) and B{p^). It is found that the extrapolations associated with 
the usage of Euclidean metric are reliable up to 1 GeV^ for calculation of masses, and up 
to around .5GeV^ for the calculation of form factors. 
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A Numerical methods 

Solutions of integral equations are performed by first discretizing the integrals 

jdqf{q)^JZwif{qi), (48) 

where Wi are integration weights for grid points g^. In order to map the grid points and 
weights from interval (—1, 1) to (0, oo) we use the arctangent mapping(Ref. [^, |3T[| ) 

^^ ^ "^^^^^1 + ^ % tan(f(l + x))' ^^'^ 



where 



It follows that 



^max J^^vned 



y{-l) = Rmin, 2/(0) = Rmed, 2/(1) = Rmax- (51) 

Therefore, one can safely control the range(i?mm, Rmax) and distribution(/2med) of grid 
points. With this discretization procedure, continuous integral equations are transformed 
into nonsingular matrix equations. 

A.l Dyson-Schwinger Equation 

The DS equation involves two unknown functions A,B which appear in two coupled 
equations(Eq. ^. After discretizing the associated integrals, one has the following matrix 
equations 

B = /iJ + G1F1, 

A = I + G2F2. (52) 

A,B,Fi, and F2 are n dimensional vectors where 

F,{{) ^ B{i)l{B{tfq} + A{t)% 
F2{i) ^ A{i)/{B{^fq^ + A{^fl 
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and Gi and G2 are n x n matrices. Coupled equations(Eq. |5^) are solved for A, B by 
forward iteration. An arbitrary initial guess for functions(vectors) A and B is entered 
on the right hand side and the resulting vectors are iteratively used for the same process 
until a stable solution is achieved. Grid points in momentum space have been chosen 
for momenta between Rmin = MeV and Rmax = lO^o" MeV where a is the relevant 
momentum scale of the problem. The median of the grid point distribution was Rmed = 5cr 
MeV. This uneven distribution of grid points ensures that the concentration of grid points 
for lower momenta, that is where the integrand is maximum, is higher. Due to the smooth 
nature of A and B functions only 40 grid points sufficed to find stable solutions. The 
number of iterations needed to find a stable result is around 20. 

A. 2 Bethe-Salpeter Equation: Inverse Iteration Method 



Here we outline the inverse iteration method originally developed in Refs. ||30| , pi 
The BS equation can be brought into the following form 

[Hm2 - «]$ = 0, (53) 

One starts with an arbitrary vector x° 

N 

X° = Ec.'^. (54) 

where $,, i = 1..N, satisfy 

[Hm2 - iu,]^, = 0, (55) 

where cUj, i = 1..N are eigenvalues of the H^^ matrix. Next, an arbitrary first guess for 

a is chosen. It should be emphasized that eigenvalues which are not equal to a have no 

physical meaning, for they do not correspond to a solution of the BS equation(Eq. |53D . 

In order to see if the initial guess a corresponds to one of the eigenvalues cji, we construct 

the operator 

K = —^ • (56) 

Um2 — a 

Operating K on state x° n times produces 

X'' = K-X' = j:^^^^^■ (57) 
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When the number of iterations n is sufficiently large(usually around ten), the dominant 
contribution to x" comes from the eigenvector $j whose eigenvalue uij satisfies |u;j — a| < 
\uji — a\ for all i = l..j — l,j + 1..N. Therefore, 



Cj 



X"+^ ^ ^—X^- (58) 

Using the eigenvector x", which is proportional to $j, one can also find eigenvalue Uj by 

u;. = ^ /^ . (59) 

X X^ 

If cuj is close enough to a, then one has a self consistent solution. This method has 
the benefit of directly singling out the eigenvalue closest to the initial guess, rather than 
finding the largest eigenvalue as in the case of straight forward iteration. There is only 
one matrix inversion involved. Distribution of the grid points in momentum space is done 
by the arctangent mapping, as in the solution of the Dyson-Schwinger equation. The 
typical number of momentum space grid points used in order to obtain stable solutions is 
around 40. For angular integrals, 20 grid points which are linearly distributed in interval 
(0, tt) proved satisfactory. 

A. 3 Confinement test 

Since the test of confinement involves a highly oscillatory integral Eq. |2B|, we have used 
a large number, 30,000, of linearly distributed grid points. The upper limit of the mo- 
mentum space integral, which is highly convergent, is 400cr. These choices allow one to 
calculate the fourier transform(Eq. ^) confidently within the timeframe shown in Fig- 
ures mi, and 0. 
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